#############################################
## PROGRAM NAME: 203d_res_muni.R           ##
## AUTHOR: MATT MLECZKO                    ##
## INPUTS:                                 ##
##         201_af_muni_fin.Rda             ##
##         202_muni_results_supp_1.qs -    ##
##         202_muni_results_supp_100.qs    ##
##                                         ##
## OUTPUTS:                                ##
##         Figures 4b-6b                   ##
##         203_muni_res1.Rda -             ##
##         203_muni_res9.Rda               ##
##         203_m1_muni_noimp.Rda -         ##
##         203_m9_muni_noimp.Rda           ##
##                                         ## 
## PURPOSE: Plot muni results              ##
##          (Figures 4b-6b)                ##
##                                         ##
## LIST OF UPDATES:                        ##
#############################################

#log <- file(# USER DEFINED PATH AND FILE NAME HERE #)
#sink(log, append=TRUE)
#sink(log, append=TRUE, type="message")

## load libraries ##
library(tidyverse) 
library(qs)
library(scales) 
library(ggtext) 


data_path <- # USER DEFINED PATH HERE #
setwd(data_path)

load("201_af_muni_fin.Rda")

## count non-missing obs ##

sum(!is.na(af.muni.fin$mgr_npov_2022))
sum(!is.na(af.muni.fin$mgr_hpov_2022))
sum(!is.na(af.muni.fin$mgr_hpov_cc_2022))

sum(!is.na(af.muni.fin$rb_npov_2022))
sum(!is.na(af.muni.fin$rb_hpov_2022))
sum(!is.na(af.muni.fin$rb_hpov_cc_2022))

sum(!is.na(af.muni.fin$mpv_npov_2022))
sum(!is.na(af.muni.fin$mpv_hpov_2022))
sum(!is.na(af.muni.fin$mpv_hpov_cc_2022))

## function to extract results ##

get_res <- function(i,j){
  
  inres <- qread(paste0(data_path,
                        "202_muni_results_supp_",
                        i,
                        ".qs"))
  
  inres.fin <- inres[[j]]
  return(inres.fin)
  
}

muni.res1 <- list()
muni.res2 <- list()
muni.res3 <- list()
muni.res4 <- list()
muni.res5 <- list()
muni.res6 <- list()
muni.res7 <- list()
muni.res8 <- list()
muni.res9 <- list()

muni.res1 <- mapply(get_res,
                    1:100,
                    1,
                    SIMPLIFY = F)

save(muni.res1, 
     file = "203_muni_res1.Rda")

muni.res2 <- mapply(get_res,
                    1:100,
                    2,
                    SIMPLIFY = F)

save(muni.res2, 
     file = "203_muni_res2.Rda")

muni.res3 <- mapply(get_res,
                    1:100,
                    3,
                    SIMPLIFY = F)

save(muni.res3, 
     file = "203_muni_res3.Rda")

muni.res4 <- mapply(get_res,
                    1:100,
                    4,
                    SIMPLIFY = F)

save(muni.res4, 
     file = "203_muni_res4.Rda")

muni.res5 <- mapply(get_res,
                    1:100,
                    5,
                    SIMPLIFY = F)

save(muni.res5, 
     file = "203_muni_res5.Rda")

muni.res6 <- mapply(get_res,
                    1:100,
                    6,
                    SIMPLIFY = F)

save(muni.res6, 
     file = "203_muni_res6.Rda")

muni.res7 <- mapply(get_res,
                    1:100,
                    7,
                    SIMPLIFY = F)

save(muni.res7, 
     file = "203_muni_res7.Rda")

muni.res8 <- mapply(get_res,
                    1:100,
                    8,
                    SIMPLIFY = F)

save(muni.res8, 
     file = "203_muni_res8.Rda")

muni.res9 <- mapply(get_res,
                    1:100,
                    9,
                    SIMPLIFY = F)

save(muni.res9, 
     file = "203_muni_res9.Rda")

## compile the results ## 

ex.msm.res <- function(i, inlist) {
  
  res <- inlist[[i]][1]
  
  return(res)
  
}

ex.noadj.res <- function(i, inlist) {
  
  res <- inlist[[i]][3]
  
  return(res)
  
}

ex.adl.res <- function(i, inlist) {
  
  res <- inlist[[i]][4]
  
  return(res)
  
}

ex.fe.res <- function(i, inlist) {
  
  res <- inlist[[i]][5]
  
  return(res)
  
}

## load results ## 

load("203_muni_res1.Rda")
load("203_muni_res2.Rda")
load("203_muni_res3.Rda")
load("203_muni_res4.Rda")
load("203_muni_res5.Rda")
load("203_muni_res6.Rda")
load("203_muni_res7.Rda")
load("203_muni_res8.Rda")
load("203_muni_res9.Rda")

##
## msm results ## 
##

m1.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.res1)

m1.msm.res.ov <- bind_rows(m1.msm.res) %>%
  filter(term == "thist")

m2.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.res2)

m2.msm.res.ov <- bind_rows(m2.msm.res) %>%
  filter(term == "thist")

m3.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.res3)

m3.msm.res.ov <- bind_rows(m3.msm.res) %>%
  filter(term == "thist")

m4.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.res4)

m4.msm.res.ov <- bind_rows(m4.msm.res) %>%
  filter(term == "thist")

m5.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.res5)

m5.msm.res.ov <- bind_rows(m5.msm.res) %>%
  filter(term == "thist")

m6.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.res6)

m6.msm.res.ov <- bind_rows(m6.msm.res) %>%
  filter(term == "thist")

m7.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.res7)

m7.msm.res.ov <- bind_rows(m7.msm.res) %>%
  filter(term == "thist")

m8.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.res8)

m8.msm.res.ov <- bind_rows(m8.msm.res) %>%
  filter(term == "thist")

m9.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.res9)

m9.msm.res.ov <- bind_rows(m9.msm.res) %>%
  filter(term == "thist")


##
## noadj results ## 
##

m1.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.res1)

m1.noadj.res.ov <- bind_rows(m1.noadj.res) %>%
  filter(term == "thist")

m2.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.res2)

m2.noadj.res.ov <- bind_rows(m2.noadj.res) %>%
  filter(term == "thist")

m3.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.res3)

m3.noadj.res.ov <- bind_rows(m3.noadj.res) %>%
  filter(term == "thist")

m4.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.res4)

m4.noadj.res.ov <- bind_rows(m4.noadj.res) %>%
  filter(term == "thist")

m5.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.res5)

m5.noadj.res.ov <- bind_rows(m5.noadj.res) %>%
  filter(term == "thist")

m6.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.res6)

m6.noadj.res.ov <- bind_rows(m6.noadj.res) %>%
  filter(term == "thist")

m7.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.res7)

m7.noadj.res.ov <- bind_rows(m7.noadj.res) %>%
  filter(term == "thist")

m8.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.res8)

m8.noadj.res.ov <- bind_rows(m8.noadj.res) %>%
  filter(term == "thist")

m9.noadj.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.res9)

m9.noadj.res.ov <- bind_rows(m9.noadj.res) %>%
  filter(term == "thist")


##
## adl results ## 
##

m1.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.res1)

m1.adl.res.ov <- bind_rows(m1.adl.res) %>%
  filter(term == "zri_2022_fin_st")

m2.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.res2)

m2.adl.res.ov <- bind_rows(m2.adl.res) %>%
  filter(term == "zri_2022_fin_st")

m3.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.res3)

m3.adl.res.ov <- bind_rows(m3.adl.res) %>%
  filter(term == "zri_2022_fin_st")

m4.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.res4)

m4.adl.res.ov <- bind_rows(m4.adl.res) %>%
  filter(term == "zri_2022_fin_st")

m5.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.res5)

m5.adl.res.ov <- bind_rows(m5.adl.res) %>%
  filter(term == "zri_2022_fin_st")

m6.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.res6)

m6.adl.res.ov <- bind_rows(m6.adl.res) %>%
  filter(term == "zri_2022_fin_st")

m7.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.res7)

m7.adl.res.ov <- bind_rows(m7.adl.res) %>%
  filter(term == "zri_2022_fin_st")

m8.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.res8)

m8.adl.res.ov <- bind_rows(m8.adl.res) %>%
  filter(term == "zri_2022_fin_st")

m9.adl.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.res9)

m9.adl.res.ov <- bind_rows(m9.adl.res) %>%
  filter(term == "zri_2022_fin_st")

## fixed effects ##

m1.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.res1)

m1.fe.res.ov <- bind_rows(m1.fe.res) %>%
  filter(term == "zri_fin")

m2.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.res2)

m2.fe.res.ov <- bind_rows(m2.fe.res) %>%
  filter(term == "zri_fin")

m3.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.res3)

m3.fe.res.ov <- bind_rows(m3.fe.res) %>%
  filter(term == "zri_fin")

m4.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.res4)

m4.fe.res.ov <- bind_rows(m4.fe.res) %>%
  filter(term == "zri_fin")

m5.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.res5)

m5.fe.res.ov <- bind_rows(m5.fe.res) %>%
  filter(term == "zri_fin")

m6.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.res6)

m6.fe.res.ov <- bind_rows(m6.fe.res) %>%
  filter(term == "zri_fin")

m7.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.res7)

m7.fe.res.ov <- bind_rows(m7.fe.res) %>%
  filter(term == "zri_fin")

m8.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.res8)

m8.fe.res.ov <- bind_rows(m8.fe.res) %>%
  filter(term == "zri_fin")

m9.fe.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.res9)

m9.fe.res.ov <- bind_rows(m9.fe.res) %>%
  filter(term == "zri_fin")

## plot the results ##

msm.res.in <- data.frame(var = c("Median gross rent (nat)",
                                 "Median gross rent (hpov)",
                                 "Median gross rent (hpov cc)",
                                 "Rent burden (nat)",
                                 "Rent burden (hpov)",
                                 "Rent burden (hpov cc)",
                                 "MPV (nat)",
                                 "MPV (hpov)",
                                 "MPV (hpov cc)"),
                         coeff = c(mean(m1.msm.res.ov$estimate),
                                   mean(m2.msm.res.ov$estimate),
                                   mean(m3.msm.res.ov$estimate),
                                   mean(m4.msm.res.ov$estimate),
                                   mean(m5.msm.res.ov$estimate),
                                   mean(m6.msm.res.ov$estimate),
                                   mean(m7.msm.res.ov$estimate),
                                   mean(m8.msm.res.ov$estimate),
                                   mean(m9.msm.res.ov$estimate)),
                         lb = c(quantile(m1.msm.res.ov$estimate,0.025),
                                quantile(m2.msm.res.ov$estimate,0.025),
                                quantile(m3.msm.res.ov$estimate,0.025),
                                quantile(m4.msm.res.ov$estimate,0.025),
                                quantile(m5.msm.res.ov$estimate,0.025),
                                quantile(m6.msm.res.ov$estimate,0.025),
                                quantile(m7.msm.res.ov$estimate,0.025),
                                quantile(m8.msm.res.ov$estimate,0.025),
                                quantile(m9.msm.res.ov$estimate,0.025)),
                         ub = c(quantile(m1.msm.res.ov$estimate,0.975),
                                quantile(m2.msm.res.ov$estimate,0.975),
                                quantile(m3.msm.res.ov$estimate,0.975),
                                quantile(m4.msm.res.ov$estimate,0.975),
                                quantile(m5.msm.res.ov$estimate,0.975),
                                quantile(m6.msm.res.ov$estimate,0.975),
                                quantile(m7.msm.res.ov$estimate,0.975),
                                quantile(m8.msm.res.ov$estimate,0.975),
                                quantile(m9.msm.res.ov$estimate,0.975)))

## no adj ## 

noadj.res.in <- data.frame(var = c("Median gross rent (nat)",
                                   "Median gross rent (hpov)",
                                   "Median gross rent (hpov cc)",
                                   "Rent burden (nat)",
                                   "Rent burden (hpov)",
                                   "Rent burden (hpov cc)",
                                   "MPV (nat)",
                                   "MPV (hpov)",
                                   "MPV (hpov cc)"),
                           coeff = c(mean(m1.noadj.res.ov$estimate),
                                     mean(m2.noadj.res.ov$estimate),
                                     mean(m3.noadj.res.ov$estimate),
                                     mean(m4.noadj.res.ov$estimate),
                                     mean(m5.noadj.res.ov$estimate),
                                     mean(m6.noadj.res.ov$estimate),
                                     mean(m7.noadj.res.ov$estimate),
                                     mean(m8.noadj.res.ov$estimate),
                                     mean(m9.noadj.res.ov$estimate)),
                           lb = c(quantile(m1.noadj.res.ov$estimate,0.025),
                                  quantile(m2.noadj.res.ov$estimate,0.025),
                                  quantile(m3.noadj.res.ov$estimate,0.025),
                                  quantile(m4.noadj.res.ov$estimate,0.025),
                                  quantile(m5.noadj.res.ov$estimate,0.025),
                                  quantile(m6.noadj.res.ov$estimate,0.025),
                                  quantile(m7.noadj.res.ov$estimate,0.025),
                                  quantile(m8.noadj.res.ov$estimate,0.025),
                                  quantile(m9.noadj.res.ov$estimate,0.025)),
                           ub = c(quantile(m1.noadj.res.ov$estimate,0.975),
                                  quantile(m2.noadj.res.ov$estimate,0.975),
                                  quantile(m3.noadj.res.ov$estimate,0.975),
                                  quantile(m4.noadj.res.ov$estimate,0.975),
                                  quantile(m5.noadj.res.ov$estimate,0.975),
                                  quantile(m6.noadj.res.ov$estimate,0.975),
                                  quantile(m7.noadj.res.ov$estimate,0.975),
                                  quantile(m8.noadj.res.ov$estimate,0.975),
                                  quantile(m9.noadj.res.ov$estimate,0.975)))

## adl ## 

adl.res.in <- data.frame(var = c("Median gross rent (nat)",
                                 "Median gross rent (hpov)",
                                 "Median gross rent (hpov cc)",
                                 "Rent burden (nat)",
                                 "Rent burden (hpov)",
                                 "Rent burden (hpov cc)",
                                 "MPV (nat)",
                                 "MPV (hpov)",
                                 "MPV (hpov cc)"),
                         coeff = c(mean(m1.adl.res.ov$estimate),
                                   mean(m2.adl.res.ov$estimate),
                                   mean(m3.adl.res.ov$estimate),
                                   mean(m4.adl.res.ov$estimate),
                                   mean(m5.adl.res.ov$estimate),
                                   mean(m6.adl.res.ov$estimate),
                                   mean(m7.adl.res.ov$estimate),
                                   mean(m8.adl.res.ov$estimate),
                                   mean(m9.adl.res.ov$estimate)),
                         lb = c(quantile(m1.adl.res.ov$estimate,0.025),
                                quantile(m2.adl.res.ov$estimate,0.025),
                                quantile(m3.adl.res.ov$estimate,0.025),
                                quantile(m4.adl.res.ov$estimate,0.025),
                                quantile(m5.adl.res.ov$estimate,0.025),
                                quantile(m6.adl.res.ov$estimate,0.025),
                                quantile(m7.adl.res.ov$estimate,0.025),
                                quantile(m8.adl.res.ov$estimate,0.025),
                                quantile(m9.adl.res.ov$estimate,0.025)),
                         ub = c(quantile(m1.adl.res.ov$estimate,0.975),
                                quantile(m2.adl.res.ov$estimate,0.975),
                                quantile(m3.adl.res.ov$estimate,0.975),
                                quantile(m4.adl.res.ov$estimate,0.975),
                                quantile(m5.adl.res.ov$estimate,0.975),
                                quantile(m6.adl.res.ov$estimate,0.975),
                                quantile(m7.adl.res.ov$estimate,0.975),
                                quantile(m8.adl.res.ov$estimate,0.975),
                                quantile(m9.adl.res.ov$estimate,0.975)))


fe.res.in <- data.frame(var = c("Median gross rent (nat)",
                                "Median gross rent (hpov)",
                                "Median gross rent (hpov cc)",
                                "Rent burden (nat)",
                                "Rent burden (hpov)",
                                "Rent burden (hpov cc)",
                                "MPV (nat)",
                                "MPV (hpov)",
                                "MPV (hpov cc)"),
                        coeff = c(mean(m1.fe.res.ov$estimate),
                                  mean(m2.fe.res.ov$estimate),
                                  mean(m3.fe.res.ov$estimate),
                                  mean(m4.fe.res.ov$estimate),
                                  mean(m5.fe.res.ov$estimate),
                                  mean(m6.fe.res.ov$estimate),
                                  mean(m7.fe.res.ov$estimate),
                                  mean(m8.fe.res.ov$estimate),
                                  mean(m9.fe.res.ov$estimate)),
                        lb = c(quantile(m1.fe.res.ov$estimate,0.025),
                               quantile(m2.fe.res.ov$estimate,0.025),
                               quantile(m3.fe.res.ov$estimate,0.025),
                               quantile(m4.fe.res.ov$estimate,0.025),
                               quantile(m5.fe.res.ov$estimate,0.025),
                               quantile(m6.fe.res.ov$estimate,0.025),
                               quantile(m7.fe.res.ov$estimate,0.025),
                               quantile(m8.fe.res.ov$estimate,0.025),
                               quantile(m9.fe.res.ov$estimate,0.025)),
                        ub = c(quantile(m1.fe.res.ov$estimate,0.975),
                               quantile(m2.fe.res.ov$estimate,0.975),
                               quantile(m3.fe.res.ov$estimate,0.975),
                               quantile(m4.fe.res.ov$estimate,0.975),
                               quantile(m5.fe.res.ov$estimate,0.975),
                               quantile(m6.fe.res.ov$estimate,0.975),
                               quantile(m7.fe.res.ov$estimate,0.975),
                               quantile(m8.fe.res.ov$estimate,0.975),
                               quantile(m9.fe.res.ov$estimate,0.975)))


## set 1: median rents ## 

msm.fin.res1 <- msm.res.in %>%
  filter(var %in% c("Median gross rent (nat)",
                    "Median gross rent (hpov)",
                    "Median gross rent (hpov cc)")) %>%
  mutate(group = case_when(var == "Median gross rent (nat)" ~ "Above Average (N=1,082)",
                           var == "Median gross rent (hpov)" ~ "High Poverty (N=233)",
                           var == "Median gross rent (hpov cc)" ~ "High-Poverty Central City (N=189)"),
         var = case_when(var == "Median gross rent (nat)" ~ "Above-average 1 (N=1,082)",
                         var == "Median gross rent (hpov)" ~ "High-poverty 1 (N=233)",
                         var == "Median gross rent (hpov cc)" ~ "High-poverty cc 1 (N=189)"),
         lb = lb,
         ub = ub,
         Model = "MSM")

noadj.fin.res1 <- noadj.res.in %>%
  filter(var %in% c("Median gross rent (nat)",
                    "Median gross rent (hpov)",
                    "Median gross rent (hpov cc)")) %>%
  mutate(group = case_when(var == "Median gross rent (nat)" ~ "Above Average (N=1,082)",
                           var == "Median gross rent (hpov)" ~ "High Poverty (N=233)",
                           var == "Median gross rent (hpov cc)" ~ "High-Poverty Central City (N=189)"),
         var = case_when(var == "Median gross rent (nat)" ~ "Above-average 2 (N=1,082)",
                         var == "Median gross rent (hpov)" ~ "High-poverty 2 (N=233)",
                         var == "Median gross rent (hpov cc)" ~ "High-poverty cc 2 (N=189)"),
         lb = lb,
         ub = ub,
         Model = "NoAdj")

adl.fin.res1 <- adl.res.in %>%
  filter(var %in% c("Median gross rent (nat)",
                    "Median gross rent (hpov)",
                    "Median gross rent (hpov cc)")) %>%
  mutate(group = case_when(var == "Median gross rent (nat)" ~ "Above Average (N=1,082)",
                           var == "Median gross rent (hpov)" ~ "High Poverty (N=233)",
                           var == "Median gross rent (hpov cc)" ~ "High-Poverty Central City (N=189)"),
         var = case_when(var == "Median gross rent (nat)" ~ "Above-average 3 (N=1,082)",
                         var == "Median gross rent (hpov)" ~ "High-poverty 3 (N=233)",
                         var == "Median gross rent (hpov cc)" ~ "High-poverty cc 3 (N=189)"),
         lb = lb,
         ub = ub,
         Model = "ADL")

fe.fin.res1 <- fe.res.in %>%
  filter(var %in% c("Median gross rent (nat)",
                    "Median gross rent (hpov)",
                    "Median gross rent (hpov cc)")) %>%
  mutate(group = case_when(var == "Median gross rent (nat)" ~ "Above Average (N=1,082)",
                           var == "Median gross rent (hpov)" ~ "High Poverty (N=233)",
                           var == "Median gross rent (hpov cc)" ~ "High-Poverty Central City (N=189)"),
         var = case_when(var == "Median gross rent (nat)" ~ "Above-average 4 (N=1,082)",
                         var == "Median gross rent (hpov)" ~ "High-poverty 4 (N=233)",
                         var == "Median gross rent (hpov cc)" ~ "High-poverty cc 4 (N=189)"),
         lb = lb,
         ub = ub,
         Model = "FE")

fin.res1 <- rbind(msm.fin.res1,
                  noadj.fin.res1,
                  adl.fin.res1,
                  fe.fin.res1)

fin.res1$Model <- factor(fin.res1$Model,
                         levels = c("MSM",
                                    "NoAdj",
                                    "ADL",
                                    "FE"))

fin.res1$group <- factor(fin.res1$group,
                         levels = c("Above Average (N=1,082)",
                                    "High Poverty (N=233)",
                                    "High-Poverty Central City (N=189)")) 


levels(fin.res1$group) <- c(
  "<b>Above Average </b>(N=1,082)",
  "<b>High Poverty </b>(N=233)",
  "<b>High-Poverty<br>Central City </b>(N=189)"
)


## thanks to Google Gemini and Microsoft Copilot for en_dash formatting help below ##

en_dash_format <- function(x) {
  x <- scales::label_comma()(x)
  gsub("-", "\u2013", as.character(x))
}

ggplot(fin.res1, aes(x=var, y=coeff, color=Model)) + 
  geom_pointrange(aes(ymin=lb, ymax=ub)) + 
  geom_hline(yintercept = 0) +
  scale_y_continuous(labels = en_dash_format) +
  theme_bw(base_size = 12) + 
  theme(
    axis.text.x = element_blank(), # Hide individual x-axis labels
    axis.text.y = element_text(size = 12, colour = "black", face = "bold"),
    axis.title.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.title.y = element_text(face = "bold"),
    panel.spacing = unit(0, "lines"), 
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text.x = element_markdown(size = 12), 
    strip.text.y = element_text(size = 12, face = "bold")
  ) + 
  facet_wrap(~group, strip.position = "bottom", scales = "free_x") +
  ggtitle("") + 
  xlab("") + 
  ylab("Difference in Monthly \nMedian Gross Rent	 ($)") + 
  scale_colour_manual(values = c("royalblue","forestgreen","red","orange"))


## set 2: rent burden ##

msm.fin.res2 <- msm.res.in %>%
  filter(var %in% c("Rent burden (nat)",
                    "Rent burden (hpov)",
                    "Rent burden (hpov cc)")) %>%
  mutate(group = case_when(var == "Rent burden (nat)" ~ "Above Average (N=1,082)",
                           var == "Rent burden (hpov)" ~ "High Poverty (N=234)",
                           var == "Rent burden (hpov cc)" ~ "High-Poverty Central City (N=189)"),
         var = case_when(var == "Rent burden (nat)" ~ "Above-average 1 (N=1,082)",
                         var == "Rent burden (hpov)" ~ "High-poverty 1 (N=234)",
                         var == "Rent burden (hpov cc)" ~ "High-poverty cc 1 (N=189)"),
         lb = lb,
         ub = ub,
         Model = "MSM")

noadj.fin.res2 <- noadj.res.in %>%
  filter(var %in% c("Rent burden (nat)",
                    "Rent burden (hpov)",
                    "Rent burden (hpov cc)")) %>%
  mutate(group = case_when(var == "Rent burden (nat)" ~ "Above Average (N=1,082)",
                           var == "Rent burden (hpov)" ~ "High Poverty (N=234)",
                           var == "Rent burden (hpov cc)" ~ "High-Poverty Central City (N=189)"),
         var = case_when(var == "Rent burden (nat)" ~ "Above-average 2 (N=1,082)",
                         var == "Rent burden (hpov)" ~ "High-poverty 2 (N=234)",
                         var == "Rent burden (hpov cc)" ~ "High-poverty cc 2 (N=189)"),
         lb = lb,
         ub = ub,
         Model = "NoAdj")

adl.fin.res2 <- adl.res.in %>%
  filter(var %in% c("Rent burden (nat)",
                    "Rent burden (hpov)",
                    "Rent burden (hpov cc)")) %>%
  mutate(group = case_when(var == "Rent burden (nat)" ~ "Above Average (N=1,082)",
                           var == "Rent burden (hpov)" ~ "High Poverty (N=234)",
                           var == "Rent burden (hpov cc)" ~ "High-Poverty Central City (N=189)"),
         var = case_when(var == "Rent burden (nat)" ~ "Above-average 3 (N=1,082)",
                         var == "Rent burden (hpov)" ~ "High-poverty 3 (N=234)",
                         var == "Rent burden (hpov cc)" ~ "High-poverty cc 3 (N=189)"),
         lb = lb,
         ub = ub,
         Model = "ADL")

fe.fin.res2 <- fe.res.in %>%
  filter(var %in% c("Rent burden (nat)",
                    "Rent burden (hpov)",
                    "Rent burden (hpov cc)")) %>%
  mutate(group = case_when(var == "Rent burden (nat)" ~ "Above Average (N=1,082)",
                           var == "Rent burden (hpov)" ~ "High Poverty (N=234)",
                           var == "Rent burden (hpov cc)" ~ "High-Poverty Central City (N=189)"),
         var = case_when(var == "Rent burden (nat)" ~ "Above-average 4 (N=1,082)",
                         var == "Rent burden (hpov)" ~ "High-poverty 4 (N=234)",
                         var == "Rent burden (hpov cc)" ~ "High-poverty cc 4 (N=189)"),
         lb = lb,
         ub = ub,
         Model = "FE")

fin.res2 <- rbind(msm.fin.res2,
                  noadj.fin.res2,
                  adl.fin.res2,
                  fe.fin.res2)

fin.res2$Model <- factor(fin.res2$Model,
                         levels = c("MSM",
                                    "NoAdj",
                                    "ADL",
                                    "FE"))


fin.res2$group <- factor(fin.res2$group,
                         levels = c("Above Average (N=1,082)",
                                    "High Poverty (N=234)",
                                    "High-Poverty Central City (N=189)")) 

levels(fin.res2$group) <- c(
  "<b>Above Average </b>(N=1,082)",
  "<b>High Poverty </b>(N=234)",
  "<b>High-Poverty<br>Central City </b>(N=189)"
)

## thanks to Microsoft Copilot for the help with this formatting function ##

en_dash_no_leading0 <- function(x, digits = 2) {
  # Base numeric formatting with consistent decimals
  lab <- scales::label_number(accuracy = 1 / (10^digits), trim = TRUE)(x)
  
  # Keep 0 as "0"
  lab[x == 0] <- "0"
  
  # Remove leading 0 for positive decimals: 0.xx -> .xx
  lab <- sub("^0\\.", ".", lab)
  # Remove leading 0 for negative decimals: -0.xx -> –.xx
  lab <- sub("^-0\\.", "\u2013.", lab)
  
  # For negative integers/other negatives: leading "-" -> en dash
  lab <- sub("^-", "\u2013", lab)
  
  lab
}

ggplot(fin.res2, aes(x=var, y=coeff, color=Model)) + 
  geom_pointrange(aes(ymin=lb, ymax=ub)) + 
  geom_hline(yintercept = 0) +
  scale_y_continuous(breaks = c(-0.01, 0, 0.01, 0.02, 0.03),
                     labels = en_dash_no_leading0) +  
  theme_bw(base_size = 12) + 
  theme(
    axis.text.x = element_blank(), # Hide individual x-axis labels
    axis.text.y = element_text(size = 10, colour = "black", face = "bold"),
    axis.title.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.title.y = element_text(face = "bold"),
    panel.spacing = unit(0, "lines"), 
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text.x = element_markdown(size = 12), # Adjust group label text
  ) + 
  facet_wrap(~group, strip.position = "bottom", scales = "free_x") +
  ggtitle("") + 
  xlab("") + 
  ylab("Difference in Proportion \nRent-Burdened Households") + 
  scale_colour_manual(values = c("royalblue","forestgreen","red","orange"))

## set 3: MPV in impoverished neighborhoods ##

msm.fin.res3 <- msm.res.in %>%
  filter(var %in% c("MPV (nat)",
                    "MPV (hpov)",
                    "MPV (hpov cc)")) %>%
  mutate(group = case_when(var == "MPV (nat)" ~ "Above Average (N=1,081)",
                           var == "MPV (hpov)" ~ "High Poverty (N=227)",
                           var == "MPV (hpov cc)" ~ "High-Poverty Central City (N=184)"),
         var = case_when(var == "MPV (nat)" ~ "Above-average 1 (N=1,081)",
                         var == "MPV (hpov)" ~ "High-poverty 1 (N=227)",
                         var == "MPV (hpov cc)" ~ "High-poverty cc 1 (N=184)"),
         lb = lb,
         ub = ub,
         Model = "MSM")

noadj.fin.res3 <- noadj.res.in %>%
  filter(var %in% c("MPV (nat)",
                    "MPV (hpov)",
                    "MPV (hpov cc)")) %>%
  mutate(group = case_when(var == "MPV (nat)" ~ "Above Average (N=1,081)",
                           var == "MPV (hpov)" ~ "High Poverty (N=227)",
                           var == "MPV (hpov cc)" ~ "High-Poverty Central City (N=184)"),
         var = case_when(var == "MPV (nat)" ~ "Above-average 2 (N=1,081)",
                         var == "MPV (hpov)" ~ "High-poverty 2 (N=227)",
                         var == "MPV (hpov cc)" ~ "High-poverty cc 2 (N=184)"),
         lb = lb,
         ub = ub,
         Model = "NoAdj")

adl.fin.res3 <- adl.res.in %>%
  filter(var %in% c("MPV (nat)",
                    "MPV (hpov)",
                    "MPV (hpov cc)")) %>%
  mutate(group = case_when(var == "MPV (nat)" ~ "Above Average (N=1,081)",
                           var == "MPV (hpov)" ~ "High Poverty (N=227)",
                           var == "MPV (hpov cc)" ~ "High-Poverty Central City (N=184)"),
         var = case_when(var == "MPV (nat)" ~ "Above-average 3 (N=1,081)",
                         var == "MPV (hpov)" ~ "High-poverty 3 (N=227)",
                         var == "MPV (hpov cc)" ~ "High-poverty cc 3 (N=184)"),
         lb = lb,
         ub = ub,
         Model = "ADL")


fe.fin.res3 <- fe.res.in %>%
  filter(var %in% c("MPV (nat)",
                    "MPV (hpov)",
                    "MPV (hpov cc)")) %>%
  mutate(group = case_when(var == "MPV (nat)" ~ "Above Average (N=1,081)",
                           var == "MPV (hpov)" ~ "High Poverty (N=227)",
                           var == "MPV (hpov cc)" ~ "High-Poverty Central City (N=184)"),
         var = case_when(var == "MPV (nat)" ~ "Above-average 4 (N=1,081)",
                         var == "MPV (hpov)" ~ "High-poverty 4 (N=227)",
                         var == "MPV (hpov cc)" ~ "High-poverty cc 4 (N=184)"),
         lb = lb,
         ub = ub,
         Model = "FE")

fin.res3 <- rbind(msm.fin.res3,
                  noadj.fin.res3,
                  adl.fin.res3,
                  fe.fin.res3)

fin.res3$Model <- factor(fin.res3$Model,
                         levels = c("MSM",
                                    "NoAdj",
                                    "ADL",
                                    "FE"))


fin.res3$group <- factor(fin.res3$group,
                         levels = c("Above Average (N=1,081)",
                                    "High Poverty (N=227)",
                                    "High-Poverty Central City (N=184)")) 

levels(fin.res3$group) <- c(
  "<b>Above Average </b>(N=1,081)",
  "<b>High Poverty </b>(N=227)",
  "<b>High-Poverty<br>Central City </b>(N=184)"
)

ggplot(fin.res3, aes(x=var, y=coeff, color=Model)) + 
  geom_pointrange(aes(ymin=lb, ymax=ub)) + 
  geom_hline(yintercept = 0) +
  scale_y_continuous(labels = en_dash_format) +
  theme_bw(base_size = 12) + 
  theme(
    axis.text.x = element_blank(), # Hide individual x-axis labels
    axis.text.y = element_text(size = 12, colour = "black", face = "bold"),
    axis.title.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.title.y = element_text(face = "bold"),
    panel.spacing = unit(0, "lines"), 
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text.x = element_markdown(size = 12), # Adjust group label text
  ) + 
  facet_wrap(~group, strip.position = "bottom", scales = "free_x") +
  ggtitle("") + 
  xlab("") + 
  ylab("Difference in Median \n Property Value ($)") + 
  scale_colour_manual(values = c("royalblue","forestgreen","red","orange"))


## END OF PROGRAM ##

#sink()
